STIN100 Datasett 3: Når får barn en tarmflora som ligner på den vi ser hos voksne?
Torgeir R. Hvidsten
Oktober 2018
1 Om øvingen
Forskerne Ekaterina Avershina og Knut Rudi her på NMBU har studert tarmfloraen til 415 mødre og deres nyfødte barn - over to år. I denne øvingen skal vi analysere disse dataene, og blant annet prøve å svare på når barn får en tarmflora som ligner på den vi voksne har.
Tarmfloraen er en betegnelse for de bakterier og andre mikroorganismer som finnes naturlig i tarmkanalen vår. Les disse artikkelene for å lære mer:
Gjennom DNA/RNA-sekvensering kan vi måle hvilke bakteriearter som finnes i en prøve og kvantifisere hvor mye det er av hver art. Les denne også, og ikke gi deg før du iallefall har lest om RNA-Seq og 16S RNAet:
Det finnes tusenvis av ulike bakteriearter i vår tarmflora, så her får vi multidimensjonale data så det holder. I denne øvingen får vi bruk for hele arsenalet av data-analyseverktøy som vi har lært i dette kurset!
1 Læringsmål
Forståelse av microbiologi-begreper:
- Forstå hvordan tarmflora-dataene er samlet inn og hva de består av.
- Forstå begreper som “tarmflora”, “bakterier” og “taxon”.
Praktiske ferdigheter i multidimensjonal dataanalyse:
- Kunne hente ut meningsfulle delmengder av et stort datasett, og gjøre analyser på disse.
- Kunne beregne og tolke avstand og korrelasjon.
- Kunne utføre og tolke resultater fra prinsipal komponent analyse og “heatmaps”.
Praktiske ferdigheter i etterrettelig forskning:
- Som på de andre øvingene skal vi knitte en R Markdown-rapport der forklarende tekst, programkode og output utgjør en fortelling om dataanalysen du har gjort.
Dette er en to ukers øving med fordypning i ett datasett. Den er normert til cirka 24 timers arbeid, forelesnings- og øvingstimer inkludert.
1.1 Krav til innlevering
Det er mye lærings- og orienteringsstoff i dette dokumentet. Det som skal inn i innleveringen, er markert med rød bakgrunn.
Tidsfrist: Mandag 12. november kl. 10:00.
Teknisk:
- En .Rmd og en .html per gruppe.
- Skriv på navnene deres under “Author”.
- .Rmd skal være kjørbar såframt nødvendige pakker er installert. Det vil si at den ikke må avhenge av filbaner som bare fins hos deg selv. Det enkleste er å legge filer du avhenger av, i samme mappe som .Rmd-fila.
- Følg en stilguide for R. (Vi er ikke pirkete, bare prøv å holde koden din oversiktlig 8-)
- Oppgi kilder for løsninger du låner.
# URL i en kommentar er nok. - Bruk R markdown til å gjøre rapporten ryddig; se
RStudio > File > Help > Cheatsheets > R Markdown Cheat Sheet. (Men ikke slit deg ut. Vi tåler fint at ting er litt klumpete og rart, bare det er begripelig.)
Vitenskapelig:
- Plott skal være selvforklarende så langt det er mulig (lettlest aksemerking, enheter på aksene hvis det passer, meningsfylte navn på factor levels), bruke grafiske virkemidler som passer til stoffet, og vise dataene på en sannferdig og dekkende måte.
- Skriv tekst om hva du skal oppnå med analysen din, hovedtrekkene i hva programkoden skal utrette, og hvordan du tolker figurer og andre analyseresultater.
1.2 Forberedelser
Les artiklene det er lenket til over og installer pakkene vi skal bruke:
tidyverse, ggplot2, gplots
library(tidyverse) # Inneholder ggplot2
library(gplots) # Heatmaps
library(readxl) # Leser inn Excel-fil2 Last inn dataene
Last in filen mor_barn_tid_tarmflora.xlsx som ligger i Canvas, og spar datarammen som data. Hvor mange rader og kolonner har data?
data <- read_excel("mor_barn_tid_tarmflora.xlsx")
dim(data)[1] 1124 1277Gjør deg kjent med dataramma. View(data) gir en dårlig opplevelse her siden tabellen er så stor, men siden dette er en tibble så skrives den veldig oversiktlig til skjerm. Se på tabellen og bruk summary-funksjonen på kolonnene 2 til 6. Beskriv hva datasettet inneholder. Hva er observasjoner, variabler og verdier?
MotherChild Timepoint Methanobrevibacter Clostridiales Dialister
Length:1124 Length:1124 Min. : 0.000 Min. : 0.000 Min. : 0.000
Class :character Class :character 1st Qu.: 0.000 1st Qu.: 0.000 1st Qu.: 0.000
Mode :character Mode :character Median : 0.000 Median : 0.000 Median : 0.000
Mean : 4.532 Mean : 1.597 Mean : 2.771
3rd Qu.: 0.000 3rd Qu.: 0.000 3rd Qu.: 0.000
Max. :344.000 Max. :434.000 Max. :784.000 Noen pekepinner:
Husk: Rader kalles observasjoner. Her er disse prøver tatt fra individer. Kolonner kalles variabler. Her er disse “bakteriearter” eller taxa (mer om det under).
Dataene inneholder observasjoner fra mor og barn (variabelen “MotherChild” forteller deg hva som er hva).
Fra mammaene er det tatt prøver under graviditeten og etter tre måneder, og fra barna er det tatt prøver etter 10 dager, tre måneder, et år og to år (variabelen “Timepoint”). Og ja, dette er bæsjeprøver.
Hvert barn har en mor (og hver mor har et barn) og denne koblingen er gitt ved at de har samme ID i variabelen “MotherChildPairID”.
De resterende variablene representer “bakteriearter” og verdiene er estimat på “antall bakterier” fra “arten” i kolonnenavnet. Jeg skriver “arten” i hermetegn, for normalt vet vi ikke hvilken art vi har funnet. Vi vet kanskje bare klassen, ordenen eller familien. Du finner en fin oversikt over taksonomien til bakterier her. Vi vil referere til disse klassifiseringene som “taxon” (litt latin: ett taxon, flere taxa). Legg merke til at flere kolonner kan henvise til samme taxon og da er det lagt på et tall til slutt for å skille de. For eksempel: “Clostridiales”, “Clostridiales.1”, “Clostridiales.2”, … Dette er altså ulike arter, men vi har bare klart å finne ut at de alle tilhører klassen Clostridiales. Kanskje er de ikke-klassifiserte arter? Moderne sekvenseringsteknologi finner ny arter i et tempo som taksonomene ikke henger helt med på.
“Timepoint” og “MotherChild” har datatype char. I plot vil vi gjerne ha tidspunktene i riktig rekkefølge i “Timepoint”, og vil vi (av grunner som blir klarere sendere) ha “mother” før “child” i “MotherChild”. Gjør disse variablene om til faktorer og bruk “levels” til å fikse rekkefølgen. Skriv ut de fem første verdiene i hver variabel.
[1] Pregnancy Pregnancy Pregnancy Pregnancy Pregnancy
Levels: Pregnancy 10days 3months 1year 2years[1] mother mother mother mother mother
Levels: mother child3 La oss telle bakterier
3.1 Hvor mange ulike taxon finnes i tarmfloraen til et individ?
Om vi ønsker å finne ut hvor mange ulike taxon som finnes i tarmfloraen til et individ så kan vi summere verdiene som er større enn null i den rader. Her har vi et stort datasett som det er vanskelig å se på. Et triks er å lage seg et lekedatasett (leke_data) som man tester koden sin på først. Under har jeg lagd et datasett med tre observasjoner fra hver gruppe av observasjoner (gruppene tilsvarer alle kombinasjoner av mor, barn og tid).
leke_data <- data[c(1:3,26,28,29,190,191,198,493:495,748:750,768,771,772), 1:6]
leke_data# A tibble: 18 x 6
MotherChildPairID MotherChild Timepoint Methanobrevibacter Clostridiales Dialister
<dbl> <fct> <fct> <dbl> <dbl> <dbl>
1 23254 mother Pregnancy 49 0 0
2 23474 mother Pregnancy 9 4 0
3 22303 mother Pregnancy 0 0 85
4 21407 mother 3months 235 0 0
5 23254 mother 3months 59 0 0
6 21831 mother 3months 0 0 4
7 22622 child 10days 0 0 0
8 22596 child 10days 0 0 0
9 22751 child 10days 0 0 0
10 20933 child 3months 0 0 0
11 21574 child 3months 0 0 0
12 23440 child 3months 2 0 52
13 21834 child 1year 0 0 1
14 22032 child 1year 0 0 0
15 22821 child 1year 0 0 0
16 23150 child 2years 1 0 0
17 22078 child 2years 0 0 0
18 22858 child 2years 0 1 0Vi kan nå bruke funksjonen rowSums, sammen med det faktum at TRUE konverteres til 1 og FALSE til 0, til å regne ut antall bakterier som er tilstedet i hvert individ (husk å ekskluder de tre første kolonnene: leke_data[,-(1:3)]).
rowSums(leke_data[,-(1:3)] > 0) [1] 1 2 1 1 1 1 0 0 0 0 0 2 1 0 0 1 0 1Sjekk at dette stemmer med leke_data.
Ok, men hva nå da? Det neste naturlige skrittet er å regne ut gjennomsnittet for mor og barn ved ulike tidspunkter for å se hvordan bakteriefloraen er forskjellig hos mor og barn og hvordan den endrer seg over tid. Nå håper jeg at dere for lengst har tenkt på dplyr-funksjonene mutate, group_by, summarise, osv!
Vi introduserte også en annen funksjon på slutten av øving 4: gather. Funksjonen mutate er fin for å lage nye variabler som er funksjoner av noen få eksisterende variabler. Her ønsker vi derimot å telle opp over alle variabler (minus de tre første), og det blir litt kronglete med mutate. gather(data, key, value, ...) samler flere kolonner i én:
dataer input-dataramma.keyogvalueer navn på to nye kolonner: Den første angir hvilken kolonne en verdi er hentet fra. Den andre angir hva verdien er....sier hvilke kolonner idatavi skal bruke. Det vil si, vanligvis angir vi (med minustegn) hvilke kolonner som ikke skal “samles”.
gather er ekstremt nyttig når vi har et datasett som inneholder mange ensartede variabler. I vårt datasett inneholder alle variablene (minus de tre første) “antall bakterier” fra ett taxon.
Bruk gather til å samle alle taxon-variablene i lekedata. Bruk pipe-operatoren %>%! Om du ikke allerede er forelsket i denne operatoren så lover jeg at det er nær forestående.
leke_data %>%
gather(Taxa, AntallBakterier, -(1:3))# A tibble: 54 x 5
MotherChildPairID MotherChild Timepoint Taxa AntallBakterier
<dbl> <fct> <fct> <chr> <dbl>
1 23254 mother Pregnancy Methanobrevibacter 49
2 23474 mother Pregnancy Methanobrevibacter 9
3 22303 mother Pregnancy Methanobrevibacter 0
4 21407 mother 3months Methanobrevibacter 235
5 23254 mother 3months Methanobrevibacter 59
6 21831 mother 3months Methanobrevibacter 0
7 22622 child 10days Methanobrevibacter 0
8 22596 child 10days Methanobrevibacter 0
9 22751 child 10days Methanobrevibacter 0
10 20933 child 3months Methanobrevibacter 0
# ... with 44 more rowsVær sikker på at du skjønner hva gather har gjort. Helt sikker.
Nå kan vi gjøre det vi vil: beregne gjenomsnittlig antall taxa som finnes i observasjons-gruppene. Lage en pipe med følgende steg (de to første har du allerede gjort):
- Ta
leke_data, - samle alle taxon-variablene i to nye variabler: “Taxa” og “Antallbakterier” (
gather), - grupper på observasjon (
group_by, grupper på “MotherChildPairID”, “MotherChild” og “Timepoint”), - summer “Antallbakterier” med verdi større en 0 og kall denne variabelen AntallTaxa (
summarise) - grupper på observasjonsgruppe (grupper “MotherChild” og “Timepoint”) (
group_by) - og til slutt ta gjennomsnitt av “AntallTaxa” og kall variabelen “GjennomsnittligAntallTaxa” (
summarise).
leke_data %>%
gather(Taxa, AntallBakterier, -(1:3)) %>%
group_by(MotherChildPairID, MotherChild, Timepoint) %>%
summarise(AntallTaxa = sum(AntallBakterier > 0)) %>%
group_by(MotherChild, Timepoint) %>%
summarise(GjennomsnittligAntallTaxa = mean(AntallTaxa))# A tibble: 6 x 3
# Groups: MotherChild [?]
MotherChild Timepoint GjennomsnittligAntallTaxa
<fct> <fct> <dbl>
1 mother Pregnancy 1.33
2 mother 3months 1
3 child 10days 0
4 child 3months 0.667
5 child 1year 0.333
6 child 2years 0.667Forsikre deg om at du forstår hva pipen over gjør: se på output fra hver nye linje i pipen og sjekk at dette stemmer med hva som ligger i leke_data. Dette hadde blitt fryktelig uoversiktlig og vanskelig om vi hadde brukt hele datasettet.
Nå er vi ferdig med å leke og er klare for å se på hele datasettet. Kjør pipen over på hele datasettet. Ser du noe mønster?
# A tibble: 6 x 3
# Groups: MotherChild [?]
MotherChild Timepoint GjennomsnittligAntallTaxa
<fct> <fct> <dbl>
1 mother Pregnancy 130.
2 mother 3months 131.
3 child 10days 22.4
4 child 3months 26.4
5 child 1year 58.6
6 child 2years 87.6Gjennomsnitt gir veldig lite informasjon om hver gruppe. Dessuten er det lettere å se mønstre i et plott. Istedenfor å se på gjennomsnittet for hver gruppe ønsker vi nå å se på de underliggende tallene i et boxplot.
Lag et boxplot over antall taxa i hver observasjons-gruppe. Hvordan tolker du resultatet nå? Hint: Dropp de to siste linjene i kodeblokka over og pipe isteden resultatet inn i ggplot. Plot boxplot med “x = Timepoint” og “y = AntallTaxa”, og bruk facet_grid til å dele opp på mor og barn. Se på eksemplene (her)[https://ggplot2.tidyverse.org/reference/facet_grid.html]. Hint: Bruk scales = "free", space = "free".
Nå er det veldig fint at vi har sortert verdiene i faktorene så alt kommer i riktig rekkefølge! Plottet over bør ha lært deg noe veldig viktig om tarmfloraen hos barn og voksene: Hvis ikke må du se på plottet en gang til. Prøv å mys.
3.2 Hvilke taxa er mest vanlige?
Over telte vi opp antall bakterier som fantes i tarmen (hadde verdi > 0) til hvert individ. Vi hadde altså et individ- (observasjon-, rad-) fokus. La oss fokusere litt på taxa nå (variablene, kolonnene).
Lag en pipe med følgende steg:
- Ta
leke_data, - samle alle taxon-variablene i to nye variabler: “Taxa” og “Antallbakterier” (
gather), - grupper på “Taxa” (
group_by), - summer “Antallbakterier” med verdi større en 0 og kall denne variabelen “AntallIndivid” (
summarise)
og test den på leke_data. Stemmer tallene?
# A tibble: 3 x 2
Taxa AntallIndivid
<chr> <int>
1 Clostridiales 2
2 Dialister 4
3 Methanobrevibacter 6Bygg ut pipen over med følgende steg:
- Bruk
ggplottil å plotte et histogram over “AntallIndivid”.
Lag et histogram med hele datasettet. Igjen vil jeg at dere ikke bare skal lage plottet men også tolke resultatet. Hint: I plottet under har jeg satt antall bins til 100. 
Hvilke taxa finnes i mer enn 250 individer? Sorter lista.
# A tibble: 97 x 2
Taxa AntallIndivid
<chr> <int>
1 Bifidobacterium 1004
2 Streptococcus 893
3 Streptococcus.2 816
4 Bifidobacterium.1 755
5 Bifidobacterium.2 717
6 Lachnospiraceae.1 676
7 Faecalibacterium 667
8 Streptococcus.3 665
9 Enterobacteriaceae.2 627
10 Bifidobacterium.4 607
# ... with 87 more rowsGoogle noen av de 10 øverste navnene. Skriv litt om tre av de: Hva er dette for en baktere (1-2 setninger) og er det kjent at de finnes i tarmen til mennesker (1 setning)?
La oss gå tilbake til histogrammet over antall individ hvert taxon finnes i. Lag et histogram for hver observasjons-gruppe. Tolk! Hvorfor er noen ruter helt tomme? Hint: Bruk pipen over som lagde ett histogram, men endre group_by og legg til facet_grid til slutt. Under satte jeg antall bins til 10. 
Ut fra histogrammene over kan det se ut som det finnes bakterier som ikke forekommer i noen observasjoner. Finn ut hvor mange disse er og fjern variablene fra data. Hva er dimensjonen på ´data´ nå?
# A tibble: 1 x 1
n
<int>
1 6[1] 1124 12714 Kvantitative data - fra telling til multivariat analyse
4.1 Normalisering
Så lang har vi bare telt forekomsten av taxa (null versus større enn null). La oss nå ta i bruk det faktum at disse dataene er kvantitative: de måler “antall bakterier” fra hvert taxon. Jeg bruker igjen hermetegn fordi disse verdiene har samme begrensinger som RNA-Seq verdiene dere lærte om i Datasett 2-øvingen. Altså: disse verdiene er basert på å telle reads og vi har ulikt antall reads fra hver prøve. Vi trenger alså å normalisere.
Jeg spurte NMBUs fremste ekspert på analyse av denne type data, Lars Snipen. Han sendte meg følgende kode. Se på den.
# The matrix read.counts must have
# - samples in the rows
# - taxa in the columns
microCLR <- function(read.counts, n.pseudo=0){
depths <- rowSums(read.counts)
RC <- read.counts + n.pseudo
RRC <- RC/rowSums(RC)
if(n.pseudo>0){
g <- apply(RRC, 1, function(x){mean(log2(x))})
X <- log2(RRC) - g
} else {
X <- RRC
}
rownames(X) <- rownames(read.counts)
colnames(X) <- colnames(read.counts)
attr(X, "depths") <- depths
return(X)
}Her har Lars altså skrevet en funksjon selv. Dette har vi ikke lært i STIN100, men som dere selv ser så er ikke dette rakett-forskning. Ved å kjøre koden over har du tilgang til funksjonen microCLR. Den tar to argument: en tabell (på samme format som vår tabell - flaks!) og en pseudo-count. Pseudo-count er et tall vi legger til alle verdier for å unngå at noen av verdiene er null, og det gjør vi fordi logaritmen til null ikke er definert. Dette heter “pseudo”-telling fordi det ikke kommer fra “ekte telling”: et read, to read, tre … Vi ser null reads, men later som vi ser noen få likevel.
Forklar hva de tre første linjene i funksjonen gjør. Om du klarer det så skal jeg forklare resten: Om du ikke angir en pseudo-count så skjer ingenting. Om du angir en pseudo-count, så blir hver verdi “x” gjort om til log2(x) - gjennomsnittet av log2-verdien til alle verdiene i den raden/observasjone. Ved å trekke fra gjennomsnittet, blir gjennomsnittet til hver rad etter normalisering lik null - dette kalles sentrering. Vi jobber ofte med log2-verdier fordi den transformasjonen sprer dataene bedre utover tallinjen og gjør verdiene mer normalfordelt (i alle fall noen ganger). Bruk funksjonen til å normalisere våre data med n.pseudo = 1, og legg resultatet i data.norm. Hint: Husk å ta bort de tre første variablene og deretter legge de tilbake etterpå (bruk bind_cols).
Lag et histogram over alle verdiene i data.norm. Diskuter hvordan resultatet stemmer med forklaringen av Lars sin funksjon som du og jeg skrev over. Hint: gather! Og under har jeg satt bins = 100. 
4.2 Prinsipal komponent analyse (PCA)
Kjør en PCA. Bruk funksjonen prcomp og legg resultatet i variablen pca.res. Lag deretter en ny tabell data.pca med prinsipal komponentene som variabler. Husk å overføre de tre første variablene i data.norm til data.pca. Dette har vi gjort flere ganger tidliger i kurset. Se for eksempel øvingen i uke 3.
pca.res <- prcomp(data.norm[,-(1:3)])
data.pca <- bind_cols(data.norm[,1:3], as.tibble(pca.res$x))Nå kan vi for alvor begynne å se på data i all sin prakt: Plott observasjonene mot de to første prinsipal komponentene. Farg punktene etter “MotherChild”. Hva ser du? 
Spennende! Det er to hypoteser om når barn får en “voksen” tarmflora: under/rett etter fødselen eller betydelig senere. La oss finne ut!
Farg punktene etter “Timepoint” og ha ulike former på punktene for mor og barn (“MotherChild”). Bruk plottet til å svare på spørsmålet i tittelen til denne øvingen: Når får barn en tarmflora som ligner på den vi ser hos voksne?.
Vi ser altså et klart mønster når vi ser på en projeksjon av hele datasettet. La oss se om vi kan finne ut hvilke taxa som er viktigst for dette mønsteret. Finn de 10 viktigste taxa for prinsipal komponent 1. Se også øvingen i uke 3 der vi diskuterte dette nærmere.
loadings <- bind_cols(Taxa=rownames(pca.res$rotation), as.tibble(pca.res$rotation))
loadings %>%
select(Taxa, PC1) %>%
arrange(desc(abs(PC1)))# A tibble: 1,268 x 2
Taxa PC1
<chr> <dbl>
1 Lachnospiraceae.1 -0.231
2 Faecalibacterium.1 -0.202
3 Faecalibacterium -0.194
4 Ruminococcaceae.3 -0.190
5 Ruminococcaceae.2 -0.171
6 Dialister.1 -0.169
7 Bacteroides.8 -0.156
8 Ruminococcaceae -0.152
9 Ruminococcus.1 -0.148
10 Bifidobacterium.8 0.147
# ... with 1,258 more rowsGoogle noen av navnene. Skriv litt om tre av de: Hva er dette for en bakterie (1-2 setninger) og er det kjent at den finnes i tarmen til mennesker (1 setning)? Kjenner du igjen noen av bakteriene fra de som fantes i flest individer (flere enn 250 individer over)? Diskuter.
Finn de 10 viktigste taxa for prinsipal komponent 2. Plott observasjonene mot det viktigste taxon i PC1 (x-akse) og det viktigste taxon i PC2 (y-akse). Hva ser du?
# A tibble: 1,268 x 2
Taxa PC2
<chr> <dbl>
1 Blautia 0.264
2 Coprococcus.3 0.261
3 Erysipelotrichaceae.7 0.215
4 Ruminococcaceae -0.192
5 Granulicatella 0.182
6 Bifidobacterium.4 0.153
7 Clostridiales.14 0.152
8 Clostridiaceae.1 0.152
9 Bifidobacterium 0.134
10 Staphylococcus -0.131
# ... with 1,258 more rows
Hmmmm. Hvor mye varians forklarer egentlig disse PCene? Skriv ut “varians forklart” for de 10 første PCene og plot deretter den komulative foredelingen av varians forklart av alle PCene. Se igjen øvingen i uke 3 for mer forklaring.
eigs <- pca.res$sdev^2
var.forklart <- eigs / sum(eigs)
var.forklart[1:10] [1] 0.15260383 0.05045699 0.03315248 0.02448756 0.02146327 0.01671577 0.01584048 0.01506326 0.01371406 0.01312990ggplot(data.frame(PC = 1:length(var.forklart), VariansForklart = cumsum(var.forklart)), aes(x = PC, y = VariansForklart)) +
geom_line()
Så de to første PCene forklarer bare 20% av variansen i data, men det er alikevel nok til å gi oss et klart bilde av hvordan tarmfloraen hos mor og barn utvikles over tid. Vi ser derimot ikke noe mønster om vi ser på enkelt-taxa, selv ikke når disse er de viktigste enkelt-taxa i PC1 og PC2. Hvordan tolker du dette?
4.3 Heatmap
Nå har jeg tenkt å slutte og holde dere i handa. Lag et heatmap. Test hva som fungerer best av korrelasjon og avstand og ulike metoder for å lage likhetstreet (hclust-funksjonen). Diskuter og tolk!. Se uke 3.
Hints:
- Husk å ta bort de tre første variablene!
- Bruk
x = as.matrix(data.norm[,-(1:3)])iheatmap.2.as.matrixer ofte svaret på feilmeldingen “must be a numeric matrix”. - Du beregner avstand/korrelasjon mellom alle par av over 1000 vektorer. Det tar et par minutter. Ikke vær så utålmodig.
- Få alt til å fungere med et lekedatasett først. Da slipper du å sitte å vente så mye:
data.norm <- data.norm[1:10,1:13].
Under ser dere mine plott, men dere behøver ikke repodusere disse eksakt. Den eneste forskjellen på de to plottene er at jeg har lagd en vertikal “colorbar” som viser hva som er mor/barn (“MotherChild”) i det første plottet og hvilke tidspunkt observasjonen kommer fra i det andre (“Timepoint”). Se argumentet RowSideColors i heatmap.2. Her er litt kode for å hjelpe dere:
# En funksjon som produserer fargene ggplot bruker. Da blir
# mor/barn og tidspunktene representert med samme farger som
# i plottene tidligere i øvingen.
gg_color_hue <- function(n) {
hues = seq(15, 375, length = n + 1)
hcl(h = hues, l = 65, c = 100)[1:n]
}
# For mor/barn:
gg_farger = gg_color_hue(2)
farger <- gg_farger[data.norm$MotherChild]
# For tidspunkt:
gg_farger = gg_color_hue(6)
farger <- gg_farger[data.norm$Timepoint]
# Du kan nå bruke "RowSideColors = farger" i heatmap.2Til sist: Jeg har også lagd en tegnforklaring til colorbar´ene. Etter at dere har kjør heatmap.2 kan det bruke legend for å legge til en tegnforklaring. Her er et eksempel:
legend(x = -0.1, y = 0.85, xpd = TRUE,
legend = levels(data.norm$MotherChild),
col = gg_farger,
pch = 15,
cex = 1,
bty = "n"
)
# OBS: x- og y-koordinatsystemet oppfører seg forskjelling i
# RMarkdown-output og i plot-vinduet i RStudio.
# Koordinatene over fungerte for meg i RMarkdown, mens i plotte-
# vinduet kommer tegnforklaringen delvis ovenpå heatmapet.
5 Avsluttende prosjekt (eksamen)
Om dere har lyst til å analysere dette datasettet videre i det avsluttende prosjektet så har jeg endel ideer for hva dere kan gjøre videre.
- Analyser mor-barn par. Variabelen “MotherChildPairID” har vi stort sett ignorert i denne øvingen. Har mor-barn en likere tarmflora en mammaer og barn generelt?
- Hvilke bakterier ligger bak de ulike klyngene av observasjoner i heatmap’et? Hent ut klynger av taxa og studer disse.
- …